Predictability of extreme events in a branching diffusion model 
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We propose a framework for studying predictability of extreme events in complex systems. Major 
conceptual elements — direct cascading or fragmentation, spatial dynamics, and external driving — 
are combined in a classical age-dependent multi-type branching diffusion process with immigration. A 
complete analytic description of the size- and space-dependent distributions of particles is derived. We 
then formulate an extreme event prediction problem and determine characteristic patterns of the system 
behavior as an extreme event approaches. In particlular, our results imply specific premonitory deviations 
from self-similarity, which have been heuristically observed in real-world and modeled complex systems. 
Our results suggest a simple universal mechanism of such premonitory patterns and natural framework for 
their analytic study. 

PACS numbers: 89.75.Hc, 89.75.-k, 91.30.pd, 02.50.-r, 91.62.Ty, 64.60.Ht 
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INTRODUCTION 

Extreme events (also called critical transitions, disasters, 
catastrophes and crises) are a most important yet least un- 
derstood feature of many natural and human-made pro- 
cesses. Among examples are destructive earthquakes, El- 
Ninos, economic depressions, stock-market crashes, and 
major terrorist acts. Extreme events are relatively rare, and 
at the same time they inflict a lion's share of the damage 
to population, economy, and environment. Accordingly, 
studying the extreme events is pivotal both for fundamen- 
tal predictive understanding of complex systems and for 
disaster preparedness (see Jill21 and references therein). 

In this paper we work within a framework that empha- 
sizes mechanisms underlying formation of extreme events. 
Prominent among such mechanisms is direct cascading or 
fragmentation. Among other applications, this mechanism 
is at the heart of the study of 3D turbulence yfl. A statis- 
tical model of direct cascade is conveniently given by the 
branching processes; they describe populations in which 
each individual can produce descendants (offsprings) ac- 
cording to some probability distribution. A branching pro- 
cess may incorporate spatial dynamics, several types of 
particles (multi-type processes), age-dependence (random 
lifetimes of particles), and immigration due to external 
driving forces |4j]. 

In many real-world systems, observations are only possi- 
ble within a specific domain of the phase space of a system. 
Accordingly, we consider here a system with an unobserv- 
able source of external driving ultimately responsible for 
extreme events. We assume that observations can only be 



made on a subspace of the phase space. The direct cas- 
cade (branching) within a system starts with injection of 
the largest particles into the source. These particles are di- 
vided into smaller and smaller ones, while spreading away 
from the source and eventually reaching the subspace of 
observations. An important observer's goal is to locate the 
external driving source. The distance between the observa- 
tion subspace and the source thus becomes a natural con- 
trol parameter. An extreme event in this system can be 
defined as emergence of a large particle in the observation 
subspace. Clearly, as the source approaches the subspace 
of observation, the total number of observed particles in- 
creases, the bigger particles become relatively more fre- 
quent, and the probability of an extreme event increases. 
In this paper, we give a complete quantitative description of 
this phenomenon for an age-dependent multi-type branch- 
ing diffusion process with immigration in W 1 . 



It turns out that our model closely reproduces the major 
premonitory patterns of extreme events observed in hierar- 
chical complex systems. Extreme events in such systems 
are preceded by transformation of size distribution in the 
permanent background activity (see e.g., [Q]]). In particu- 
lar, general activity increases, in favor of relatively strong 
although sub-extreme events. That was established first by 
analysis of multiple fracturing and seismicity H Si, and 
later generalized to socio-economic processes 13]. Our re- 
sults suggest a simple universal mechanism of such pre- 
monitory patterns. 
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MODEL 

The system consists of particles indexed by their gener- 
ation k = 0, 1, Particles of zero generation (immi- 
grants) are injected into the system by an external forcing. 
Particles of any generation k > are produced as a result 
of splitting of particles of generation k — 1. Immigrants 
(k = 0) are born at the origin x := (xi, . . . , x n ) = ac- 
cording to a homogeneous Poisson process with intensity 
(i. Each particle lives for some random time r and then 
transforms (splits) into a random number (3 of particles of 
the next generation. The probability laws of the lifetime r 
and branching j3 are rank-, time-, and space-independent. 
New particles are born at the location of their parent at the 
moment of splitting. 

The lifetime distribution is exponential: P{r < t} = 
1 — e _At , A > 0. The conditional probability that a par- 
ticle transforms into n > new particles (0 means that it 
disappears) given that the transformation took place is de- 
noted by p n . The probability generating function for the 
number (3 of new particles is thus 

h( S )=J2PnS n . CD 

n 

The expected number of offsprings (also called the branch- 
ing number) is B : = E(f3) = h'(\) (see e.g., 0]). 

Each particle diffuses in W l independently of other par- 
ticles. This means that the density p(x, y, t) of a particle 
that was born at instant at point y solves the equation 

X p = DA xP (2) 



dt \ ^ dx? 



with the initial condition p(x, y , 0) 
lution of (0 is given by 



<5(x — y). The so- 



p(x,y,t) = (AirDt) " /2 exp 



x 



(3) 



ADt 

where |x| 2 = ^\ xf. 

It is convenient to introduce particle rank r : = r max — k 
for an arbitrary integer r max and thus consider particles of 
ranks r < r max . This reflects our focus on direct cas- 
cading, which often assumes that particles with larger size 
(e.g., length, volume, mass, energy, momentum, etc.) split 
into smaller ones according to an appropriate conservation 
law. Figure [T]illustrates the model dynamics. 



SPATIO-TEMPORAL PARTICLE RANK 
DISTRIBUTIONS 

The model decsribed above is a superposition of inde- 
pendent branching processes generated by individual im- 
migrants. We consider first the case of a single immigrant; 
then we expand these results to the case of multiple immi- 
grants. Finally, we analyze the rank distribution of parti- 
cles. Proofs of all statements will be published in a forth- 
coming paper. 
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D) Disiance lo the origin |x|=0 
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FIG. 1: Example of a 3D model population. Different panels 
show 2D subspaces of the model 3D space at different distances 
|x| to the origin. Model parameters are fi = X = 1, D = 1, 
B = 2. Circle size is proportional to the particle rank. Different 
shades correspond to populations from different immigrants, the 
descendants of earlier immigrants have lighter shade. The clus- 
tering of particles is explained by the splitting histories. Note 
that, as the origin approaches, the particle activity significantly 
changes, indicating the increased probability of an extreme event. 



Single immigrant 

Let Pk,i(G,y,t) be the conditional probability that at 
time t > there exist i > particles of generation k > 
within spatial region G C W given that at time a sin- 
gle immigrant was injected at point y. The corresponding 
generating function is 



F k (G,y,t;s) = Y / PkAG,y 1 t)s\ 



(4) 



Proposition 1 The generating functions F^(G, y, t; s) 
solve the following recursive system of non-linear partial 
differential equations: 

J^-Ffc = -DA y F k — \F k + \h (F k _i) , k > 1,(5) 

with the initial conditions F k (G, y, 0; s) = 1, k > 1, and 
F (G,y,t;s) = (l-P)+Ps, (6) 
where P := e~ At J G p(x, y, t)dx.. 

Next, consider the expected number A k (G, y , t) of gen- 
eration k particles at instant t within the region G produced 
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by a single immigrant injected at point y at time t = 0. It 
is given by the following partial derivative (see e.g., ||4j]) 



(7) 



Consider also the expectation density A k (x, y, t) that sat- 
isfies, for any GcK", 

A k (G,y,t) = f A fc (x,y,i)dx. (8) 

JG 

Corollary 2 The expectation densities A k (~x,y,t) solve 
the following recursive system of linear partial differential 
equations: 

dA k 



dt 



D A x A k - XA k + XBA k ^, k > 1, (9) 



with the initial conditions A k (x., y, 0) = 0, k > 1, 

A(x,y,0) = <5(y-x), 

A)(x,y,i) = e~ A V(x,y,i), t > 0. (10) 

The solution to this system is given by 

A k (x,y,t) = (A f ^ A)(x,y,t). (11) 



The system (|9]) has a transparent intuitive meaning. The 
rate of change of the expectation density A k (x, y, t) is af- 
fected by the three processes: diffusion of the existing par- 
ticles of generation k in M n (first term in the rhs of Q), 
splitting of the existing particles of generation k at the rate 
A (second term), and splitting of generation k — 1 particles 
that produce on average B new particles of generation k 
(third term). 



Multiple immigrants 

Here we expand the results of the previous section to the 
case of multiple immigrants that appear at the origin ac- 
cording to a homogeneous Poisson process with intensity 
u. The expectation A k of the number of particles of gener- 
ation k is given, according to the properties of expectations, 
by 

A(x,t) = f A k (x,0,s)uds (12) 
Jo 

The steady-state spatial distribution A k (x) corresponds to 
the limit t — > oo and is given by 

Here z := |x| ^f\/D, v = k — n/2 + 1 and K v is the 
modified Bessel function of the second kind. 



Rank distribution and spatial deviations 

Recall that the particle rank is defined as r = r max — k. 
The spatially averaged steady-state rank distribution is a 
pure exponential law with index B: 



A, 



A k (x,0,t)ndtdx = ^B k oc B~ r . (14) 
A 



To analyze deviations from the pure exponent, we con- 
sider the ratio 7fc(x) between the number of particles of 
two consecutive generations: 

A(x) 



7fc(x) := 



(15) 



A+i(x) ' 

For the purely exponential rank distribution, ^4fc(x) = 
cB k , the value of 7fc(x) = 1/B is independent of k and 
x; while deviations from the pure exponent will cause 7 fc 
to vary as a function of k and/or x. Combining (PT3l and 
(TT5T ) we find 

2(fc + l) K v (z) 

7fc(x) - 



Bz K„ +1 {z) 



(16) 



where, as before, z := |x| \J\jD andv = k- n/2 + 1. 

Proposition 3 The asymptotic behavior of the function 
7^ (z) is given by 

( oo, v < 0, 

lim 7fc (z) = <^ l_ f l+ n_\ „ (17) 



7k(z) 



B V 2u 
2{k + 1 
Bz 



, u>0, 
oo, fixed k, 



(18) 



1 / Tl \ 

^ (l + ^j , fc^oo, fixed .2.(19) 



Proposition [3] allows one to describe all deviations of 
the particle rank distribution from the pure exponential law 
(fj~4b - Figure |2 illustrates our findings. First, Eq. (fl~9T > im- 
plies that at any spatial point, the distribution asymptoti- 
cally approaches the exponential form as rank r decreases 
(generation k increases). Thus the deviations can only 
be observed at the largest ranks (small generation num- 
bers). Analysis of the large-rank distribution is done using 
Eqs. ([T7ll,(fl"8l. Near the origin, where the immigrants enter 
the system, Eq. ( fTTl l implies thatj k (z) > j k+ i(z) > 1/B 
for v > 0. Hence, one observes the upward deviations 
from the pure exponent: for the same number of rank 
r particles, the number of rank r + 1 particles is larger 
than predicted by (fT4l) . The same behavior is in fact ob- 
served for v < (the details will be published else- 
where). In addition, for v < the ratios ^ k {z) do not 
merely deviate from 1 / B, but diverge to infinity at the ori- 
gin. Away from the origin, according to Eq. ( fl8t . we have 
7fe(z) < 7fc+i(£) < 1/B, which implies downward devia- 
tions from the pure exponent: for the same number of rank 
r particles, the number of rank r + 1 particles is smaller 
than predicted by ( fT4l . 
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FIG. 2: Deviations from self-similarity: Expected number A k {z) 
of generation k particles at distance z from the origin (cf. Propo- 
sition [3). The distance z is increasing (from top to bottom line 
in each panel) as z — 10~ 3 , 2, 5, 10, 20. Model dimension is 
n = 1 (panel A), n = 3 (panel B), n = 5 (panel C), and n — 10 
(panel D). Other model parameters: /x = X = 1, D = 1, B = 2, 
r max = 21. One can clearly see the transition from downward 
to upward deviation of the rank distributions from the pure expo- 
nential form as we approach the origin. 



DISCUSSION 



shots of a 3D model at different distances from the source. 
One can see that, as the source approaches, the following 
changes in the background activity emerge: a) The inten- 
sity (total number of particles) increases; b) Particles of 
larger size become relatively more numerous; c) Particle 
clustering becomes more prominent; d) The correlation ra- 
dius increases; e) Coherent structures emerge. In other 
words, the model exhibits a broad set of premonitory phe- 
nomena previously observed heuristically in real and mod- 
eled systems: multiple fracturing [6], seismicity [5], socio- 
economics [7], percolation @], hydrodynamics, hierarchi- 
cal models of extreme event development [1]. These phe- 
nomena are at the heart of earthquake prediction algorithms 
well validated during 20 years of forward world-wide tests 
(see e.g., H). 

In this paper we analyse only the first-moment proper- 
ties of the system; such properties can explain the premon- 
itory intensity increase (item a above) and transformation 
of the particle rank distribution (item b). At the same time, 
the framework developed here allows one to quantitatively 
analyze other premonitory phenomena; this can be readily 
done by considering the higher-moment properties. 

This study was supported in part by NSF grant No. ATM 
0620838. 



Motivation for this work is the problem of prediction of 
extreme events in complex systems. Our point of depar- 
ture is a classical model of spatially distributed population 
of particles of different ranks governed by direct cascade 
of branching and external driving. In the probability the- 
ory this model is known as the age-dependent multi-type 
branching diffusion process with immigration |4J]. We in- 
troduce here a new approach to the study of this process. 
We assume that observations are only possible on a sub- 
space of the system phase space while the source of ex- 
ternal driving remains unobservable. The natural question 
under this approach is the dependence of size-distributions 
of particles on the distance to the source. The complete 
analytical solution to this problem is given by the Proposi- 
tion Q] 

It is natural to consider rank as a logarithmic measure 
of the particle size. If we assume a size-conservation law 
in the model, the exponential rank distriburtion derived in 
( TBI) corresponds to a self-similar, power-law distribution 
of particle sizes, characteristic for many complex systems. 
Thus, the Proposition [3] describes space-dependent devia- 
tions from the self-similarity (see also Fig. 12); in particular, 
deviations premonitory to an extreme event. The numerical 
experiments (that will be published elsewhere) confirm the 
validity of our analytical results and asymptotics in a finite 
model. 

The model studied here exhibits very rich and intriguing 
premonitory behavior. Figure Q] shows several 2D snap- 
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